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f**^ ^ Abstract. In the past several decades there have been a number of proposals for computing with dual 

, forms of non-abelian Yang-Mills theories on the lattice. Motivated by the gauge-invariant, geometric 

picture offered by dual models and successful applications of duality in the U(l) case, we revisit the 
question of whether it is practical to perform numerical computation using non-abelian dual models. 
Specifically, we consider three-dimensional SU{2) pure Yang-Mills as an accessible yet non-trivial 
case in which the gauge group is non-abelian. Using methods developed recently in the context of 
spin foam quantum gravity, we derive an algorithm for efficiently computing the dual amplitude and 
describe Metropolis moves for sampling the dual ensemble. We relate our algorithms to prior work in 
non-abelian dual computations of Hari Dass and his collaborators, addressing several problems that 
have been left open. We report results of spin expectation value computations over a range of lattice 
sizes and couplings that are in agreement with our conventional lattice computations. We conclude 
with an outlook on further development of dual methods and their application to problems of current 
interest. 
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1. Introduction 



cn 

^ I In this work we describe recent progress in lattice gauge theory (LGT) computations in a model 

' dual to pure Yang-Mills theory. The dual model we consider is based upon the character expansion of 

\ the amplitude at each fundamental plaquette, a procedure that has been known since the early days of 

■ lattice gauge theory [25, 36, 42]. To date, the most common use of this duality transformation has been 
in the strong-coupling expansion, as described for example in [22, 26, 32] and references therein. Because 

' of this historical association with strong-coupling approximations, we should emphasize here that it is 

I the exact dual model that we compute with, sampling from the full space of dual configurations to find 

' the expectation values of observables. The first steps towards such dual computations in the non-abelian 

case were taken by Hari Dass et al. in [16, 17, 20]; more recently, interest in dual computations of lattice 
. , ■ Yang-Mills theory has grown within the spin foam community [10, 33, 39]. 

r> ' Perhaps the most significant change in passing to the dual model is that the variables include discrete 

■ labels assigned to lattice plaquettes. This aspect alone gives a very different character to the computation 
when compared to the conventional formulation of LGT, where continuous group-valued variables are 
assigned to lattice edges. Another significant difference is the fact that the discrete labels must satisfy 
certain constraints, making the choice of Metropolis moves non-trivial (see Section 2). The allowed 
configurations can be viewed as closed, branched surfaces colored by the irreducible representations of 
the gauge group. To illustrate, some renderings of non-zero dual configurations are shown in Figure 1. 

As general motivation for the dual approach, we observe that the dual configurations have a gauge- 
independent meaning rooted in the statistical geometry of the two-dimensional branched surfaces, pro- 
viding a compelling geometric picture for the evolution of the physical degrees of freedom. For example, 
in the strong-coupling limit, the dual picture provides a straightforward analytic proof of confinement 
(see, e.g., [32] and references therein). By making topological excitations manifest, a dual model is well 
suited to evaluating proposed mechanisms of quark confinement, such as the dual superconductor pic- 
ture. This advantage has already been demonstrated in the abelian C/(l) case [23, 37, 38, 40, 44, 45, 46] 
and has also been remarked upon in the non-abelian case [33]. Recently, effective theories have been 
derived directly from the dual theory [11], and so we believe a computational framework for numerical 
computations within the dual theory is timely. 



While results have been reported with dual models in the U{1) theory, the non-abelian case has 
presented a greater computational challenge. As a testing ground for dual non-abelian simulation 
algorithms, we focus in the present work on SU (2) pure Yang-Mills in three space-time dimensions. 
The observable we study is the average spin (see Section 3) which is convenient for testing as it takes a 
simple form on both sides of the duality. In describing our algorithm and results, we address potential 
problems for non-abelian dual computations that were discussed in the work of Hari Dass et al. [17, 20], 
the most critical of which were the construction of ergodic moves and a sign problem. The results 
obtained with our Metropolis algorithm appear to overcome these difficulties for the range of lattice 
sizes and coupling constants presented here. With the current algorithm, a difficulty (discussed in 
Section 4.2) emerges at weak coupling, however our expectation is that it will be resolved with a more 
refined form of the algorithm. 

Although there is an intriguing possibility that dual algorithms may eventually outperform their con- 
ventional equivalents in certain contexts, the present work focuses on the general features of computing 
with dual algorithms and verifying their correctness against a conventional lattice code. In future work, 
we will address more optimized implementations of the algorithms described here, with the ultimate 
goal of providing results in four dimensions with lattice sizes typical of contemporary lattice QCD. An 
important next step will be to extend our methods to Wilson loop observables crucial to the study of 
quark confinement and glueball spectra; this work is currently in progress [7] . 

The paper is organized as follows. In Section 2, we define the dual model under study and describe 
a Metropolis algorithm for numerical computations. In Section 3, we describe computations on a 2^ 
lattice that allow us to verify threefold agreement (within statistical error) between dual Metropolis, 
conventional Metropolis and the exact partition function. In Section 4, we describe our computations 
for lattice sizes of 4'^, 8'^, and 16^ and report agreement between dual and conventional Metropolis 
results, again within statistical error. Section 5 provides a brief outlook on generalizing the algorithm 
and Section 6 presents our conclusions. Some details relevant for performing the computations are given 
in the Appendix, which also shows how the dual amplitude formula used by Hari Dass [17, 20], originally 
given in [3] and discussed in [13, 14], fits into the spin foam formalism used to develop our algorithm. 




Figure 1. Two examples of dual configurations on a 6'^ lattice. Opposite sides of the 
lattice are identified. In the right-hand figure a sheet of flux spanning the lattice is 
present, allowed due to the non-trivial homology of the 3-torus. In the left-hand figure 
there are no such sheets present. 

2. Review of The Dual Model And Algorithm Construction 

In this section we define the dual model for SU{2) pure Yang-Mills on a hyper-cubic lattice and 
describe a Metropolis algorithm for performing dual computations. Of particular relevance to our 
present work is the application of diagrammatic methods and spin foam formalism, as described for 
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example in [10, 33, 34], a critical feature being the expression of group integrals of representation 
matrices and their contractions as sums over intcrtwincrs^. While the spin foam formalism for dual 
models is a more recent development, the construction of exact dual non-abelian models has a longer 
history [1, 2, 3, 14, 15, 18, 19]. Traditional derivations of the duality transformation and their use in 
strong-coupling expansions can be found in the texts [22, 32]. 

2.1. Review of pure Yang-Mills theory on the lattice. First we recall the Euclidean partition 
function of pure Yang-Mills theory in D dimensions, with gauge group G = SU{N), where N > 2 (we 
shall later specialize to the SU{2) case). It takes the form 

(1) ^ " / '^^ ^^P("'^)' 

with ^J' the gauge field, S the action functional, and VA the functional integration measure. In the 
continuum version of the theory the standard action functional is 

(2) s^S[A] = ^ J d^xF;:^Fr, 

where is the field strength tensor and g the continuum coupling. Unfortunately, the continuum 
functional measure T>A is not well-defined. 

One way to give the above path integral rigorous meaning and, at the same time, make it amenable 
to computational treatment, is to put the theory on a discrete finite lattice. The simplest variant uses a 
hyper-cubic lattice. Let E and P denote respectively the sets of edges and plaquettcs of a hyper-cubic 
lattice in D dimensions. The gauge field A is replaced by gauge group elements ge assigned to each 
oriented lattice edge e E E. The same edge with opposite orientation gets g~^ instead of g^- The 
functional integral measure can now be replaced by an integral over the product of |£^| copies of G using 
Haar measure: 

(3) = n "^de. 

At the same time, the action functional is replaced by a discretized version, S = S[g], that must 
reproduce the continuum action S[A] as the lattice spacing is taken to zero. The discretized action is 
usually split into a sum over plaquettcs, S[g] = X^pep ^i9p)j where the group element gp is the holonomy 
around an oriented plaquette p. That is, gp = gig^g^gAi where gi is either the group element assigned 
to the ith edge of p or its inverse if the orientations of p and the ith edge are opposing. This yields the 
conventional lattice partition function 

(4) Z = I JJrfffee-^feP^(S''\ 

There are many candidate discretized plaquette actions S{gp). While the Wilson action [43] is perhaps 
the most well-known in conventional LGT (it was also used in the dual computations of [16, 17, 20]), 
a variety of actions S{gp) leading to the correct continuum limit are known and have been used in the 
literature [21, 28, 29]. In the present work, we use the heat kernel action [30]; in the dual model this 
action leads to plaquette factors that are particularly easy to compute. The heat kernel action (at lattice 
coupling 7) for a fundamental plaquette p and plaquette holonomy gp is 

where the heat kernel K, which is a function of a group element g and of a "time" parameter t, satisfies 
a diffusion type differential equation 

(6) ^^K{g,t) = AK{g,t), K{g,0) = Sj{g). 

Here A is the Laplace-Beltrami operator on G and dj is the delta function at the group identity /. 
The denominator in (5) represents a normalization of the partition function in which flat holonomies 



An intertwiner is a map between representations of a group tliat commutes with the action of the group. 
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{qp = I) are assigned an amplitude of unity. We shall follow the common practice of discussing the 
phase structure of a lattice theory using the /3 parameter P — 

We now turn to the definition of the dual model for the specific case of G = SU{2) pure Yang-Mills 
in three dimensions. Starting from the conventional formulation of the lattice partition function Z 
given in (4) above, the duality transformation can be applied (see Appendix A.l) to yield the following 
expression for Z in terms of the dual variables: 

j \ i vev eeE / \peP J 

Here V denotes the vertex set of the lattice, while the summations over i and j range over all possible 
edge and plaquette labellings, respectively. A plaquette labelling j assigns an irreducible representation 
of SU (2) to each element of P. These representations are labelled by non-negative half-integers (we will 
denote this set by ^N) and arc referred to as spins; a labelling j is thus a map j: P —f ^N. An edge 
labelling i, on the other hand, is valued in a basis of maps that intertwine the representations of the 
plaquettes incident on the same edge. In our present case, the choice of basis corresponds to a grouping 
of the four incident plaquette spins into two pairs. When such an edge splitting has been made, the 
intertwiners may also be labelled by spins, as described just before Definition 2.2 and in Appendix A. 2. 
Different choices of splitting can be made, but, some are more computationally efficient than others. 
In writing (7), we assume a fixed choice of splitting has been made and so an edge labelling is a map 
i: E^\n. 

In the first pair of parentheses of (7), there is a product of 18j symbols, each of which is a function of 
the 18 spins which label the 12 plaquettes and 6 edges incident to a vertex v. we denote the spins which 
appear by jy and iy . Next to it is a product of edge normalizations iV^ depending on the edge spin «e 
and on the four spins je labelling the plaquettes incident on e. It is important to recognize that the ISj 
symbol and the normalization factors N*^ are purely representation-theoretic quantities (independent 
of the action chosen) and that, from a computational viewpoint, they represent the non-trivial part of 
the amplitude evaluation. Efficient algorithms can be found (using diagrammatic techniques similar to 
those used in [8]) for computing the 18j symbols and edge normalizations. Two of these are reviewed 
in Appendix A. 2. 

In the second parentheses of (7) there is a product of factors depending on jp only; these arise from 
the character expansion coefficients of the heat kernel action (5) and are clearly straightforward to 
compute. 

For the purposes of this paper, we define a spin foam to be an assignment of spins and intertwiners 

to the plaquettes and edges of the lattice^, respectively. We define supported spin foams to be those 
with non- vanishing amplitude, and denote the set of them by JF+. In terms of the supported spin foams 
we can write (7) as: 

(8) z=Y^ A{f)= E (ni8r(i.,>)n^'(^-j^)"M (n«"^"'^"'^'n2ip+i)) . 

feJ^+ .feT+ \vev eeE I \veP j 

It turns out that the spin foams that actually make a non-zero contribution to Z are highly constrained, 
a fact that is not manifest in (7) where every possible edge and plaquette coloring contributes a term. 
Moreover, once the plaquette spins have been specified, there is a limited range of intertwiner spins that 
give a non-zero amplitude. The details of these constraints will be discussed in Section 2.2. 

The essential challenge for performing Metropolis simulations with this model is to find a set of moves 
that connect all the supported spin foams . To see why this is a non-trivial task, observe that if one 
has any supported spin foam and adds or removes a half-unit of spin from a single plaquette, the result 
will be a spin foam with vanishing amplitude; this follows immediately from the parity condition (see 
Definition 2.1 below), enforced by the 18j symbol. Changing single plaquettes by two half-units of spin 
preserves parity but restricts one to a single parity class. Thus a move involving a single plaquette is 



In spin foam quantum gravity, the 2-cells and 1-ceIIs of a simplicial complex are typically used, rather than the 
plaquettes and edges of a cubic lattice. The general definition of spin foam was introduced by Baez in [4] . 
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"too local" and so moves involving multiple plaquettes and edges simultaneously are needed. We now 
turn to this problem. 

2.2. Algorithm Construction — Ergodic Moves. We first define the set of admissible plaquette 
configurations. These configurations are precisely those plaquette labellings that satisfy edge admissi- 
bility at every edge of the lattice: 

Definition 2.1 (Edge admissibility). The spins assigned to plaquettes incident to an edge are said to 
be edge admissible if the parity and triangle inequality conditions are satisfied. Writing ji, j2, js and 
ji for the four spins incident to a given edge, these conditions are: 

(1) Parity: 

ji + J2 + J3 + J4 is an integer. 

(2) Triangle Inequality: for each permutation {k, I, m, n) of {1, 2, 3, 4} we have 

3k+3l +3m > 3n- 

These conditions are equivalent to the existence of a non-zero invariant vector in the SU{2) representa- 
tion ji oi2 oi3 oi4- 

If a splitting has been made for the edge, say with ji and j-^ on one side and js and ji on the other, 
then the triangle inequality can be written in the less symmetric form 

(9) |ji - J2 1 < is + J4, lis -ji\< h + J2- 

Representation-theoretically, the less symmetric form says that there exists a non-trivial intertwiner 
between ji ® and ja ® ji. A basis of such intertwiners is labelled by spins which match parity with 
ji +j2 and j3 + j4 and fall into the specific range defined below. This leads us to the following definition 
of an admissible spin foam: 

Definition 2.2 (Admissible spin foam). A spin foam is admissible if and only if: 

(1) The assignment of spins to plaquettes is everywhere edge admissible. 

(2) For every edge e G E, i^ satisfies 

ie + ji + 32 and ie + js + ji are integers 

and 

ie&Ie= [|il - j2\ ,il +i2] H [|j3 - ji\ ,^3+^4], 

where [a, b] = {j € \ a < j < b} . 

We call the interval Ie the range of the triangle inequalities at edge e. It is guaranteed to be non-empty 
by edge admissibility. We denote the set of admissible spin foams by T"^. 

At this point we observe that while any supported spin foam is admissible, the converse does not 
necessarily hold, as the amplitude can vanish despite achnissibility conditions being satisfied. We de- 
fine exceptional spin foams to be admissible spin foams of vanishing amplitude and denote them by 
= . In practice, exceptional spin foams are rarely encountered during simulation. If one 

assumes that exceptional spin foams are sufficiently isolated (i.e. do not form surfaces that separate the 
configuration space) then the question of ergodicity will not be affected, i.e. it will be sufficient to show 
ergodicity of moves on admissible spin foams. In what follows, we shall make this assumption and refer 
to it as the isolation of exceptionals hypothesis. 

Before considering ergodicity with respect to spin foams, it is usehil to first describe moves by which 
any admissible plaquette configuration can be reached from any other admissible plaquette configuration: 

Definition 2.3 (Plaquette cube move). A plaquette cube move consists of 

(1) A choice of 3-cell of the cubic lattice. These are simply the smallest cubes of the lattice, with 6 
plaquettes, 12 edges, and 8 vertices. 

(2) For each of the 6 plaquettes of the 3-cell, a change is made to that spin by -t-5 or — | units of 
spin. 
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It can be shown that these moves are ergodic on the space of admissible plaquette configurations 
(ignoring for the moment non- trivial global topology of the lattice^); a similar argument is used for 
instance in [5]. 

We next observe that, within our present context where the amplitude depends on intertwiner la- 
bellings, the plaquette cube moves will transform between admissible plaquette eonfigurations but inad- 
missible spin foams, unless some simultaneous changes are made to the intertwiners labelling the edges 
of the cube. This leads us to define spin foam cube moves as follows: 

Definition 2.4 (Spin foam cube move). A spin foam cube move consists of the following: 

(1) A plaquette cube move. 

(2) For each edge in the 3-cell chosen for the plaquette cube move, a randomly chosen change is 
made to the intertwiner labelling that edge, from a fixed set of possible changes. The possible 
changes are such that for any plaquette cube move, there is always a change that results in an 
admissible spin foam. The intertwiner changes that are proposed depend on how the edges are 
split within the selected cube; further details can be found in Appendix A. 3. 

We should mention here that for all the spin foam moves discussed in this section there will be many 
moves that are rejected immediately due to violation of the constraints; these are nonetheless counted 
as moves in order to satisfy detailed balance. 

By construction, spin foam cube moves will reach some subset of intertwiner labellings for each 
plaquette configuration. The possibility that this may omit some admissible intertwiner configurations 
leads us to introduce an independent edge move. 

Definition 2.5 (Spin foam edge move). A spin foam edge move consists of the following: 

(1) An edge is selected. 

(2) The intertwiner label is incremented or decremented by two half-units of spin (to preserve 
parity). 

To see that the combination of spin foam cube and edge moves are ergodic on admissible spin foams 
(up to non-trivial global homology), we argue as follows. Observe that the spin foam cube moves allow 
us to move freely amongst the admissible plaquette eonfigurations. To reach any given spin foam, one 
first obtains the associated plaquette configuration through the spin foam cube moves. We then apply 
whatever spin foam edge moves are needed to set the edge labels to their given values. This is possible 
because the admissible ranges le are connected by spin foam edge moves. 

We now turn to the the issue of configurations arising from the non-trivial global topology of the 
lattice. It is common in lattice computations to impose periodic boundary conditions which are toroidal; 
i.e. in the present case the lattice is a discretization of the 3-torus. For these bomidary conditions, if 
one introduces a sheet of half- unit spin having non- trivial global topology (see Figure 1, right), the spin 
foam cube moves are not able to remove it. Such a sheet (and its deformations by cube moves) can be 
introduced in any of the three directions of the 3-torus, in correspondence with the three generators of 
the second homology group of the 3-torus. A simple way to move between configurations corresponding 
to different homology classes is to create and remove sheets of half-unit spin; we do this by introducing 
the following moves: 

Definition 2.6 (Spin foam homology move). A spin foam homology move consists of the following: 

(1) A 2-dimensional plane of plaquettes spanning the lattice is selected. Applying periodic boundary 
conditions, this plane topologically defines a 2-torus wrapped around the 3-torus. 

(2) For the selected plane, change each edge and plaquette label contained in it by the same amount, 
either +^ or — i. 

We note that the cost of computing the change in amplitude induced by these moves scales with 
(where L is the side length of the lattice), making them costly for large lattices; we shall return to this 
point in our discussion below. 

Now we have all the components to give a simple statement of our algorithm. 

■^If boundary conditions result in a lattice of non-trivial homology, homology changing moves must be added, as 
described below. 
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Algorithm 2.7. (Ergodic spin foam algorithm) A single iteration of the algorithm consists of the 
following: 

(1) Apply a single spin foam cube move, spin foam edge move, or spin foam homology move. 

(2) For the edges affected by the move, immediately reject the move if any spin foam admissibility 
condition is violated. 

(3) Accept the move if > a, where A^ew and A^id denote the spin foam amplitudes before 
and after the move, and a is a sample drawn from from the uniform distribution on the unit 
interval. 

Step 1 gives a choice from among a set of moves that are ergodic on admissible spin foams, not just within 
a homology class but between classes as well. For each type of move, detailed balance is guaranteed 
as any move and its inverse is proposed with equal probability. The relative probability of proposing 
a cube, edge, or homology move can be freely chosen, for example to tune the algorithm for better 
convergence or acceptance rate. Step 2 could in principle be absorbed into Step 3, as the zeroes of the 
18j symbols enforce all the admissibility constraints. In practice, admissibility is very easy to check and 
so one can often avoid computing the full 18j symbol. Step 3 is just the usual Metropolis acceptance 
condition; obscrvablcs depending on the dual variables can be measured at each iteration or at regular 
intervals. The absolute value of the amplitudes in Step 3 are necessary due to the possibility of negative 
amplitudes (arising from the 18 j symbols). The standard sign trick for computing expectation values 
in this situation is reviewed in Section 2.3.2. 

In summary, we have defined a dual Metropolis algorithm for non-abelian SU{2) pure Yang-Mills 
in the case of three dimensions. Assuming the isolation of exceptionals hypothesis, this algorithm is 
ergodic on supported spin foams. Before presenting results for our algorithm in Sections 3 and 4, we 
review some practical points of implementation. 

2.3. Algorithm Construction — Practical considerations. When using the above algorithm, 
there are a number of novel considerations that don't apply to the conventional case. We discuss 
some of the more important ones here. 

2.3.1. Spin foam versus plaquette configuration simulation. It is possible to define an alternative dual al- 
gorithm based upon plaquette configurations (rather than spin foams, which include intertwiner labels), 
using the plaquette cube moves to move ergodically between all admissible plaquette configurations. By 
applying the intertwiner expansion and fast evaluation of the 18j symbols, one can automate the com- 
putation of plaquette configuration amplitudes. This would proceed by writing the partition function 
(7) in terms of an intertwiner independent amplitude: 

(10) z = ^i(j), i(j) ^ ( E n ^^f^iv,3v) n N-{ie,3e)-'] I n ^'^'"'^'"^^x^jp +m, 

j \ i veV eGB / \peP J 

where the amplitude A(j) is a function of plaquette labelling j only; i.e. the full sum over intertwiner 
labels has been absorbed into the amplitude for a given plaquette labelling. 

The question then becomes whether this algorithm is practical. We found that the answer depends 
crucially on the strength of the coupling. For comparison, we implemented the plaquette configuration 
simulation and were able to achieve good agreement with conventional results at sufficiently strong 
coupling (/? < 1.8), where configurations involving large intertwiner STims are penalized by the amplitude. 
However, starting at /3 w 1.8 the simulation began to spend time in plaquette configurations with a large 
space of admissible intertwiner labellings; because the summations are nested there are a huge number 
of edge labelled terms being evaluated for a single plaquette configiuation. Compounding this situation 
is a second critical drawback, in that for many configurations (and increasingly so at weaker coupling) 
the terms in the intertwiner summations couple edges that are arbitrarily distant in the lattice. Thus, 
at large lattice sizes and weak coupling the cost of computing the change in amplitude due to a single 
move will on average be of the order of the lattice volume. 

The spin foam algorithm defined in Section 2.2 avoids both problems by making local changes to both 
plaquettes and edge labels. In summary, rather than importance sampling the plaquette configurations 
with a generally expensive, non-local amplitude (the intertwiner sum) , one applies importance sampling 
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to the innermost terms of (7), which carry both plaquctte and edge labels. While this involves a larger 
space of configurations, the evaluation of amplitudes proceeds extremely rapidly and independent of 
coupling. The observation that such a local algorithm exists for dual Yang-Mills theories (and also for 
dimensions higher than 3) was made by Halliday and Suranyi in [14], although they do not make explicit 
use of the spin foam formalism. 

2.3.2. Sign Problem. Because the 18j symbol can take on negative values, the most straightforward 
form of the Metropolis algorithm is not applicable. Nonetheless, one can try to avoid this problem by 
applying what is commonly referred to as the "sign trick", e.g. [41]. We note that the sign trick was 
also employed in the work of Hari Dass [16]. Letting e(/) denote the sign of the amplitude A{f) for 
spin foam /, one can express the expectation value of an observable 0{f) as follows: 

/Ef=»-+ 0{f)e{f)\AU)\\ 

^^^^ _ EfeF^ 0{f)A{f) _ 0(/)e(/) \A{f)\ _ [ '%^,,\au)\ ) _ (eO),, 



\ E/<=F+I^(/)I / 
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Observe that we now have a ratio of two expectation values, which can be estimated using simulations 
of a system governed by the absolute value \A{f)\ of the original amplitude (we indicate this by adding 
a subscript ||). A common failing of this approach is that for many systems studied in quantum Monte 
Carlo (where the sign problem arises frequently), both numerator and denominator of the expression 
decay exponentially with a parameter under study. In such cases there arc more sophisticated ap- 
proaches available; invariably, resolving a serious sign problem requires exploiting foreknowledge of the 
distribution of positive and negative weights in the configuration space. 

As discussed in Section 4, for the range oi (3 < 2.85 and lattice sizes up to the maximum studied 
(16^), the numerator and denominator of (11) are of order unity. In this regime, the sign problem is 
very tame and does not prevent accurate results from being obtained. 

2.3.3. Cost of Hom,ology Moves. Due to the extended nature of the structures added and removed by 
the homology moves (involving on the order of L'^ vertices and plaquettes), the associated change in 
amplitude becomes increasingly expensive to calculate for larger lattices. The situation for the present 
case of SU(2) Yang-Mills in three dimensions can be summarized as follows. Homology changing moves 
are necessary to achieve agreement with conventional code on very small lattices, such as the 2^ lattice 
discussed in Section 3. However, we have found in practice that as the lattice size increases, the frequency 
with which homology moves are attempted can be greatly decreased and ultimately set to zero without 
affecting the result. This is because adding a sheet of half-unit spin has a very low acceptance rate for 
larger lattices. We should mention however that for different spin foam models, particularly those which 
are in a deconfining phase, the non-trivial homology sectors of the ensemble will likely make significant 
contributions for arbitrarily large lattice sizes, and so cannot be neglected. 

3. Verification of Dual Metropolis Code on A Small Lattice 

Before proceeding to simulations at larger lattice size, an implementation of the dual algorithm 
described in Section 2 was verified on a 2^ lattice. For this lattice and all lattices discussed in the paper, 
toroidal boundary conditions were applied. The observable measured was ^p^pjp^ the average spin 
per plaquette. The expectation value for this observable is: 

Af) 

Prom the translation invariance of our lattice, this is a good estimator for the expectation value of spin 
on a single plaquette. 

While too small to be of physical interest, the 2^ lattice size makes it easy to check agreement between 
different Monte Carlo algorithms, as convergence is fast on smaller lattices. Moreover, by applying a 
spin cutoff it is possible to exactly compute the (cut-off) partition function of the dual model on a 
2'^ lattice. This is done for cutoffs jcut = \ and jcut = 1, which correspond to the restrictions that 
i G {O;^} and j G {O, respectively. While jcut — 1 gives on the order of 10^° configurations 

on the 2^ lattice, the case of jcut = | already requires on the order of 10^^ configurations and is not 




practical. In Figure 2, we show the results for spin cutoff jcut = 1- For the conventional computations, 
the cut-off is realized by truncating the character expansion of the heat kernel; e.g. for jcut = 1 only the 
first three terms of the character expansion are used. This truncation leads to a mild sign problem for 
the conventional calculation that can be resolved using the sign trick for the range of f3 presented here. 
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Figure 3. Conventional and dual results for the 8^ lattice. 



The main plots in Figure 2 each show two curves which appear on top of each other. The inset plots 
show the relative errors 

(13) 0-)e -(i)c e,^%;i^, 

where the subscripts c, d and e refer to the conventional, dual and exact calculations. The expectation 
values (j)^ and {j)^ and their standard deviations were found from the results of 12 runs of 7.5 x 10^ 
simulation steps each. The error bars are 3a in magnitude. The data is consistent with the expected 
error being zero. 

4. Testing The Dual Algorithm 

4.1. Results. In this section we describe results obtained on an 8^ lattice over a range of coupling 
constants. For dual spin foam simulations, 25 runs of 10^° moves were used to generate the data shown 
in Figure 3; for the conventional results, 15 nms of 10^ moves were used. The observable measured 
was the expectation value (j) of the spin observable, as described in Section 3. In the main part 
of Figure 3, the expectation value for both conventional and dual results are shown, with standard 
deviation computed from the results of different runs. Qualitatively, wc sec that the two algorithms are 
in very good agreement, as the curves fall nearly on top of one another. To better exhibit the detailed 
behavior, the inset plots show the relative error Ccd = ^^hjy^- '^^^ standard deviation in the relative 
error was computed by combining the standard deviations of the expectation values. All error bars are 
3(7 in magnitude. We find agreement between the two algorithms within about 1%, which is within the 
estimated error. 

The results for lattices of size 4"^ and IG'' are very similar. The above computations were done with 
the ladder recoupling of the 18j symbol (see Appendix A. 2.1) due to its simpler implementation. For 
the range of coupling presented, we found that the results were insensitive to increases in the cut-off 
past jcut = 3, which was the specific value used for both the dual and conventional data. 

4.2. Slow-down at weak coupling. With the current implementation of the algorithm, there is a 
problem in obtaining results for (3 > 2.85. For these (3, the Metropolis simulation eventually becomes 
"trapped" in a region of higher spin than untrapped configuations. These trapped regions appear to 

10 



have steep or narrow exits, as the simulation remains trapped for the duration of the computation. This 
behavior is also observed in some runs at stronger coupling (small /?) ; the trapped mode generally takes 
a longer time to develop as coupling is increased, and is absent entirely from most of the runs at the 
strongest couplings considered. 

While the precise nature of this slow-down is not yet fully understood, it shows some of the symptoms 
of similar problems noted for dual U (1) in three dimensions. In [12], long autocorrelation times (not seen 
with the conventional code at the same /3) are observed past a certain /3. Slow-down of dual algorithms 
due to overly local moves at weak coupling is also alluded to in [44, 45, 46] for C/(l) in four dimensions. 

By analogy with the existing literature on the dual U{1) case, it is possible that the solution will 
involve adding moves that simultaneously change spins over regions that are larger than the minimal 
cube. A useful analog to consider may be spin systems, where the variables are also discrete and 
slow-down has been successfully addressed using cluster algorithms [27]. In our case, the presence of 
admissibility and parity constraints tends to make it difficult to find moves with reasonable acceptance 
rates. As is common with cluster algorithms, the "size" of the extended moves will likely need to be 
tuned with /3. 

5. Outlook for further Development 

5.1. Generalizations. Although we have so far limited ourselves to SU{2) pure Yang-Mills in three 
dimensions, the algorithm we describe is in principle straightforward to extend to more physically 
relevant models. 

5.1.1. Wilson Loop Observables. The theory behind Wilson loop observables in the dual spin foam 
picture can be found in the work of [33]. The basic challenge for computations is the fact that, in 
the dual, the Wilson loop is a ratio of two partition functions with different admissibility conditions. 
In a forthcoming paper, we describe how to extend the method described here to define a Metropolis 
algorithm for Wilson loop observables [7]. 

5.1.2. Four Dimensions. In four dimensions on a hyper-cubic lattice, the 18j symbol is replaced by its 
higher-dimensional analog^, a 32 j symbol (involving 24 plaquettes and 8 edges); for the heat kernel action 
the plaquette weights take the same form as in (7). As the edges will have six incident plaquettes, the 
triangle inequalities will be more involved. Assuming these differences can be accounted for, we expect 
an algorithm with the same general form will be feasible. 

Of particular interest in D = 4 will be the behavior of autocorrelation times and simulation slow- 
down as the critical point (not present in three dimensions) is approached in the continuum limit of 
large lattice and weak coupling. 

5.1.3. Higher N Lie groups and their q- deformations. The algorithm has direct generalizations to SU (3) 
and even SU{N). Vertex amplitudes (generalizations of the I8j symbol) will likely be more intricate, 
as for higher N more than one discrete variable is needed to label a given irreducible representation. 
The basic structure of the algorithm in terms of moves on plaquette and intertwiner labels should go 
through, although admissibility constraints can be expected to take a more complicated form. 

An interesting feature of the dual model is the ease with which g-deformed gauge groups SUq{N) can 
be treated once the infrastructure for computing with SU (N) is in place. One illustration of this is recent 
work [9] where computations are performed with a g-deformcd version of the Riemannian Barrett-Crane 
model of spin foam quantum gravity. In four dimensions, where one expects to approach a critical point 
in the continuum limit, it would be interesting to consider how computations of correlation functions 
are affected by g-deformation, for example whether or not distinct universality classes emerge. 



In the same way that the dual model [2] and [3] is equivalent to a particular recoupling of the 18j symbol (see 
Appendix A.2.2), we expect the 32j to agree with past work on dual forms of Yang-Mills in four dimensions [14]. 
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5.1.4. Other Extensions. The above mentioned extensions of the method are straightforward. There are 
also intriguing (and more involved) possibilities of extending the method into the realm of dynamical 
fermions and to gravity. For incorporating dynamical fermions into a dual model, we are investigating 
an approach in which the polymer expansion of the fermion determinant plays a key role. For computing 
pure Yang-Mills coupled to spin foam quantum gravity, work by Oriti and Pfciffcr [35] provides a useful 
starting point. Developments in the theory of dual models [34, 39] indicate that the dual viewpoint has 
broad applications. We hope that concrete computational results will encourage its further development 
and offer an alternative means for analyzing theories of current interest. 



6. Conclusions 

Previous work with computing in the dual of non-abelian theories due to Hari Dass et al. [16, 17, 20] 
brought to light a number of important issues, but it had not been clear how these might be fully 
resolved. 

The work reported here establishes that one can indeed obtain agreement with conventional lattice 

gauge theory using simulation with the dual form of the theory. Currently, the algorithm we describe 
has provided results on a 16"^ lattice for (3 up to 2.85, and work is underway to extend this to larger 
lattices and weaker coupling (higher (3). We expect addressing higher (3 will involve introducing moves 
that change spin foams over a larger region of the lattice, as our algorithm appears to exhibit the 
symptoms of crossover to a delocalized, disordered phase as observed in [12, 44, 45, 46] for the case of 
dual three-dimensional U{1). In four dimensions, the expectation is that similar behavior will be seen 
as the critical point is approached. The fact that cluster algorithms have been remarkably successful 
at addressing critical slow-down in discrete spin systems is promising in this regard and it is likely that 
some insights from this area may apply. 

The ultimate impact of the sign problem is still not fully understood from the present work, as the 
sign expectation begins decreasing away from unity at close to the same /? where the current simulation is 
obstructed by impractically long autocorrelation times. Assuming a practical resolution to the trapping 
problem can be found, it will be interesting to see whether or not the sign problem becomes a serious 
obstacle at weaker coupling. 

An important property of our algorithm is the use of both plaquette labels and edge (intertwiner) 
labels as variables of the simulation. This allows the Metropolis algorithm to proceed through local 
moves whose changes in amplitude can be evaluated very rapidly. Using the spin foam formalism, 
ergodicity is fairly straightforward to show; it is essentially an extension of an algorithm used in previous 
(quantum gravity) spin foam computations [5, 9] to the case where intertwiner labels are changed as 
well. Our claim of ergodicity does however depend on the unproven hypothesis that exceptional zeros 
of the 18j symbol are isolated, i.e. the zeros do not form surfaces that separate the space of admissible 
spin foams. The agreement of our results with conventional computations can be viewed as evidence 
for this hypothesis. 

Another interesting aspect of our algorithm is the necessity (for full ergodicity) of introducing moves 
that change homology class, in the case where boundary conditions lead to non-trivial global topology 
of the lattice. While for the specific model studied here such configurations have negligible weight for 
large lattices, in more general spin foam models (e.g. those in a deconfining phase) these configurations 
may play a more important role. 

Also worth emphasizing is the important role of diagrammatic techniques in generating efficient 
algorithms for the ISj symbols. In particular, the formula of Anishetty et al. [2, 3] used in [16, 17, 20] 
can be derived very economically from a particular splitting of the 18j symbol, as we show in the 
Appendix. 

While still in its very early stages (particularly in comparison with several decades development of 
conventional LGT codes), we believe dual computations have a promising future ahead. The algorithm 
we describe here immediately generalizes to higher N gauge groups and to four dimensions, and work is 
underway to incorporate dynamic fermions into the dual framework. The next major step is the com- 
putation of Wilson loop observables, currently in progress [7]. Finally, the application of the technique 
to spin-foam quantum gravity coupled to matter provides a long-term motivation for the present work. 
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Appendix A. The dual model, 18 j symbol algorithms and compatible intertwiner moves 

A.l. Derivation of the dual model. This section sketches some of the steps of the transformation 
from the conventional to the dual form of the lattice Yang-Mills partition function, (4) and (7) re- 
spectively. Our approach is inspired by the spin foam picture, and is closest to that found in [10]. 
Non-abelian dual models have also been analyzed from a spin foam perspective in [33, 34]. 

We begin by observing that, due to gauge invariance, the plaquette action S{gp) of (4) depends 
only on the conjugacy class of its argument. Thus, its exponential can be expanded in terms of group 
characters Xj 

(14) e-^(^)=^c,-x,(5), 

j 

where j ranges over the equivalence classes of finite-dimensional irreducible unitary representations of 
the gauge group G. Substituting into (4) and interchanging the order of summation and integration 
yields 

(15) I U^Sellcj^xM- 

bp} ee-E pep 

At this point it is convenient to specialize to a D = 3 cubic lattice with periodic boundary conditions 
and to fix an orientation for the plaquettes and edges of the lattice. Choose a right-handed set of xyz 
axes for the lattice. Orient all of the edges in the positive coordinate directions. Every lattice cube is 
in the first octant of one of its vertices. Take each of the three plaquettes of the cube that are incident 
to this vertex and orient it in the counterclockwise direction, as seen from outside the cube. It is easy 
to sec that this choice of orientations is translation invariant, that the orientation of each edge agrees 
with two of the four plaquettes incident on it and is opposite to the other two, and that every plaquette 
has two edges whose orientations agree with its own and two that do not. 

With this choice of orientation, the holonomy around a plaquette p is gp = gig2g3^ 9^^ , where 
51)52,53 and (?4 are the group elements associated to the edges of the plaquette p, starting with an 
appropriate edge and going cyclically. Recall that the inverse g^^ is used if the orientation of edge i 
does not agree with that of p. Thus 

(16) xjM = uMl uMh UjM'fc UjM'Td , 

where Uj{g)^ denotes a matrix element with respect to a basis of the j representation. If we insert (16) 
into (15) and collect together factors depending on the group element ge, we get a product of independent 
integrals over the group, each of the form 




Here and below we use a graphical notation for tensor contractions, defined as follows. Each wire 
represents a matrix element of the unitary representation labelling it. Parallel wires represent products 
of such matrix elements. The four matrix elements in (17) come from the characters associated to the 
four plaquettes incident on the edge e. The free ends of the wires represent the indices of these matrix 
elements. The wires can be joined together into loops, one for each plaquette. The joining corresponds 
to contracting with other matrix elements from different edge integrals to form the product of characters 
as in (16). 
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Equation (17) defines a projection operator on the space of linear maps j4 (g) js j\ ® j2- It is 
the usual group averaging operator whose image is precisely the intertwiners. Since it is a projection 
operator, it can be resolved into a sum over a basis of intertwiners /, : j4 (S) jz ji ji, 




where the intertwiners /* : ji (g) j2 — > ji (g Js are chosen such that the trace of the composite 

is zero whenever i' ^ i and non-zero if i' = i. The projection property is readily verified. 
If, for each edge of the lattice, we fix a term i in the above summation, we can contract the intertwiners 
Ii and /* with those coming from the other edges. At each vertex of the lattice, there will be six such 
intertwiners, and their contraction can be graphically represented as an octahedral network that we call 
the 18j symbol: 



+z 



(19) 




-z 



The vertices are labelled by the directions of the associated lattice edges emanating from the given 
lattice vertex, namely ±x, ±j/, and ±z. The value of the ISj symbol depends on the choice of basis 
elements Ii and /*, in (18), the six summation indices i labelling the edges, and the 12 incident plaquette 
labels j. Each normalization factor N = (I*, Ii) depends on the choice of basis elements at an edge, the 
summation index i on that edge, and the four plaquettes incident on that edge. Note that the choice of 
basis can be made independently at each edge. 

The discussion up to this point has been quite general, assuming a 3-dimensional cubic lattice. Next, 
we specialize to G = SU{2) and give the plaquette action character expansion coefficients. For the heat 
kernel (5), the expansion coefficients take the particularly simple form [30] 

(20) e-«('') = — ^{2j + l)e-^^(^+i)x, (5), j = 0, i, 1, . . . , 

K[l, 2 ) j 

where K{g, t) is defined by (6). Putting these pieces together, we obtain the dual formula for the lattice 
Yang-Mills partition function 

(21) Z = ^(5]ni8j'"(^-^'')n^'(^-j'^)"M i]l(2j, + l)e-'^i^^^''+A , 

j \ i vev eeE / \peP J 

2 

where an overall numerical factor of K{I, ^) per plaquette has been discarded. This precisely reproduces 
Equation (7), where we described the notation we arc using for the plaquette and edge labellings j and i. 

A. 2. Efficient algorithms for the 18j symbol via recoupling. In order to perform computations 
with (21), we first need to choose explicit basis elements Ii and /* of the spaces of intertwiners that 
appear in (18). Below, we consider two patterns for choosing such bases for each edge of the lattice, 
one we call the ladder recoupling and one we call the tetrahedral recoupling. They lead to different 18j 
symbols and have different properties with respect to lattice translations. 
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A. 2.1. The ladder recoupling. Recall that for compatible spins j, k and m, there is an intertwiner 
j ®k ^ m that is unique up to scale. To be explicit, we choose the specific intertwiner defined in [6, 
2.5.4], and we denote it by 



(22) 

m m 

Similarly, we use the same reference^ to define 

m 

(-) X 

j k j k 

It is well-known that for fixed ji, j2, js and ji, the intertwiners 

(24) 77 ^ ]j 

ii h 

form a basis of the space of intertwiners ji (g) js — > j\ (g) j2, as i varies over admissible spins. We call this 
the vertical splitting. 

There is also a second vertical splitting, given by interchanging and ji, which changes the inter- 
twiner by a factor of (— l)-^3+j4-* 'pjjg geometry of the lattice provides a natural way to choose between 
the two: we make sure that the plaquette labels on the left {j^ and ji above) are part of the same lattice 
cube, and same for the labels on the right. 

A convenient choice of dual basis is given by 

(25) /r = 

U ] 3 

One can check that = for i' ^ i. We next need to evaluate the normalization factor 

(26) N^ = {ir,ID 

In order to accomplish this, we now explain how to relate our tensor contraction diagrams to spin 
networks. While it would be possible to work entirely with tensor contraction diagrams, there are two 
reasons to switch to the spin network notation. First, spin networks do not require that the edges be 
directed, which relieves us of some complicated bookkeeping. Second, by using spin networks, we can 
take advantage of many existing formulas and software libraries for computing spin network evaluations. 

Recall that a spin network is a trivalent undirected ribbon graph whose edges are labelled by spins. 
One assigns a value to a spin network in the following way. First, draw it in the plane, in general 
position, with the ribbon flat. Then, read it from top to bottom, interpreting the vertices as the 
trivalent intertwiners discussed above, and interpreting cups and caps as certain intertwiners which 




'Note that our diagrams are read downwards, while those of [6] are read upwards. 
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can introduce signs. If the spin network is closed, the resulting intertwiner is a map from the trivial 
representation to itself, and so can be identified with a complex number. The result is independent 
of the embedding in the plane, a fact that is quite useful in computations. In particular, the trivalent 
intertwiners are chosen carefully so that the evaluation remains unchanged when inputs are deformed 
into outputs and vice versa. We refer the reader to [6] for more details. While the starting point is 
different, the formulas given in [24], with A = 1, also apply to these diagrams^. Note that some other 
authors have slightly different conventions, e.g. some take A = —1. 

We will now work out how to compute tensor contractions using spin networks. Take a tensor 
contraction diagram involving just the trivalent intertwiners discussed above and draw it in the plane 
such that all edges are pointing downwards except for some edges which leave the bottom of the diagram 
and loop around to reenter at the top. If we erase the arrows, the resulting spin network will have the 
same interpretation as the tensor contraction diagram, except for the signs introduced in the cups and 
caps. One can show that the difference is exactly a factor of (—1)^'^, where J is the sum of the spins 
labelling the edges that loop around. 

As a first example, the value of a loop labelled with j in the tensor notation is the dimension 2j + 1 of 
the representation. However, in the spin network notation, the value of a loop is Aj = (— l)^^ (2j + 1). 

Similarly, the value of the edge normalization factor is 



(27) 



AT" = (If 




(-1)^ 



J 2 



J 3 



iAh,j2,i)0{j3,j4,i) 



where 6{a, b, c) stands for the value of the following theta network: 

a 

(28) e{a, h,c)i ^ 



Its value is given explicitly in [24. Chapter 9]. 

Note that the conversion sign factor (— 1)'^* from (27) can be expressed as (— l)20i+i2) by appealing to 
the parity constraints. Since each plaquette is "outgoing" from two edges, each plaquette spin contributes 
(—I)'*-' = 1. In other words, the conversion sign factors from the edge normalizations iV" cancel. It can 
be shown that the conversion sign factor for the 18j symbol appearing in (19) is independent of the 
edge splitting and can be written as (— 1)^-^, where, for instance, J = j+x+y +J+x-2 + J-y-z- Each 
plaquette spin shows up in exactly one such sign factor, so the signs combine to give (—1)^'^'°', where 
Jtot is the sum of all plaquette spins. Note that on a lattice with two or more odd side-lengths, this 
sign factor can be non-trivial. 

Next we must work out the value of the ISj symbols that arise using the vertical splitting. The 
corresponding spin network is obtained by applying this splitting to the vertices of the octahedron 
shown in (19) and erasing the arrows from its edges. A method for evaluating this spin network is 
shown in Figure 4. The calculation is similar to that of [8], where a "ladder" structure also appears. 
The recoupling move 



(29) 



m 

n 



6(a, d, m,)9{b, c, m) 



m 



'Note that while we use half-integer spins, [24] uses twice-spins, which are always integers. 
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Vertex splitting 
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-x-z Redrawing in plane 



Redrawing as a ladder (wrapped in a circle) 
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Recoupling to a bubble chain 



E 



+x+y 



Hx 3+x+v rn 

i-z j-y-z j+x- 



J+v+z i+z m 

3—x—y ^—x 3—X+ 



3+x+v i+v rn 

3-x-z i-z 3+y-. 



Hz 3+x+z rn 

r-y 3-x-y 3-y+z 



»+» 3+y+z m 
I'—x 3—x—z 3— x+y 



3+x+z t+x 



-y-z ^-y 3+x-y 



n: 



••)«(•••) 

Figure 4. The ladder recoupling of the 18j symbol; a single sum. 



is applied to each of the six "rungs" of the ladder, producing a chain of bubbles. The function of six 
spin labels appearing in (29) is the tetrahedral network, shown in the last step of Figure 5. The value of 
the tetrahedral network is given explicitly in [24, Chapter 9] and is closely related to the Wigner-Racah 
6j symbol of angular momentum theory [31, Appendix B], see (31). 

Because of Schur's Lemma, the six independent sums from the recoupling moves become a single 
sum. The bubbles are proportional to the identity, weighted by a theta network divided by a loop. 
Six theta networks arising from the bubbles cancel against six of the twelve theta networks from the 
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recoupling moves to give the six theta networks shown in the final line. The bubbles also contribute six 
loop factors (A^) in the denominator, which exactly cancel the loop factors from the recoupling. The 
final result can be written as: 



(30) E 



n 



6 bubbles 



m 



9{m, i+a:,j-y-^)e{m,i-a;,j+y+z)0{m, i+y,j-^-z)e{m,i-y,j+a;+z)0{m,i+z,j-x-y)O{m,i-z,j+x+y)' 



where the arguments of the six tetrahedral networks are those that appear in the last line of Figure 4. 
The explicit relation between the tetrahedral network and the Wigner-Racah 6j symbols is 

Jl J2 J3 



(31) 



Jl 



'2 

32 
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VWuJ2,h)9{juj2,h)e{Jl,j2, J3)e{juJ2, J3)| 



Jl J2 J3 
J2 J3 



Note the row swap and the fact that the four theta networks correspond to the four triples of spins from 
the 6j's arguments that must satisfy triangle inequalities. For reference, \9{a, b, c)| = (— l)°'+''+^6'(a, b, c). 

The 18j symbol described in this section was used in computing the data appearing in Sections 3 
and 4. 

A. 2. 2. The tetrahedral recoupling. Next we consider a different splitting of the vertices of the octahedron, 
which we call the tetrahedral recoupling. The 18j symbol that arises here is more efficient to compute 
than the I8j symbol for the ladder recoupling, because it does not require a sum. However, the splitting 
is not translation invariant, which makes it slightly harder to work with. This section is not needed in 
the rest of the paper, but is useful as a comparison to other sources and will be important for future 
calculations. 

We begin by considering a different basis for the space of intertwiners j4 (8) js ji (8) j2 ■ It is given 
by the horizontal splitting 



(32) 



J4 



Jl 



J3 



as i varies over admissible spins. Note that it makes no difference which way the arrow on the edge 
labelled by i points. 

There is also a second horizontal splitting, given by interchanging J3 and J4. As we did for the vertical 
splitting, we choose between the two by requiring that ji and j\ label plaquettes that are part of the 
same cube. Unlike the vertical splitting, the two horizontal splittings are not in general related by a 
sign. 

A convenient dual basis is given by 



Jl 



(33) 

The traces , if) work out to be 

(34) N^ = {l'^,\lf) 



J4 



J2 



J3 



(-1) 



2(ii+j2) 



For the tetrahedral recoupling, we use the vertical splitting (24) on three edges and the horizontal 
splitting on the opposite edges. The pattern we use is indicated in the first step of Figure 5. The result 
is a planar network consisting of four triangles connected to one another in a tetrahedral pattern. 

For this recoupling, a minor complication arises because the intertwiner splittings are diff'erent on 
opposite edges. This means that a simple translation of the given tetrahedral 18j symbol to neighboring 

18 




-x-y 




Figure 5. The tetrahedral recoupling of the 18j symbol; no sum. 



vertices does not correspond to a consistent choice of basis for the intertwiners. This is easily overcome 
by dividing the lattice into a checkerboard of odd and even sites, and alternately nsing the original and 
reflected versions of the 18j symbol. For periodic boundary conditions, this docs limit one to lattices 
with even side-lengths, but this constraint is not serious in practice. 

As was the case with the ladder rcconpling, the conversion sign factors from the normalization factors 
cancel, and the conversion sign factors from the 18j symbols give a factor of (—1)^'^'°'. In this case, 
because the checkerboard pattern forces even side-lengths, one can show that (— l)2-^tot = i 

The diagrammatic relation 



h 




can be used to collapse each of the four triangles into 3-valent vertices, as shown in the last step of 
Figure 5. The final result is simply: 
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Figure 6. (a) The two types of splitting for an edge in a cube (side view), (b) Two 
possible splittings associated with the ladder recoupling for the edges of a cube. 



(36) 
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Because this formula is essentially a product of tetrahedral networks with no auxiliary summations, it 
is highly efficient to compute. 

The tetrahedral recoupling is easy to express in terms of Wigncr-Racah 6j symbols using relation (31). 
In order to compare our work to other work, we give some of the details. Ignoring signs for the moment, 
the theta networks from the edge normalizations (27) and (34), the conversion formula (31), and the 
vertex amplitude (36) all cancel. Now we collect the signs of the theta networks. The theta networks 
from the edge normalizations contribute a sign of (— l)2«(_i)Ji+J2+i3+j4^ where i labels the edge and the 
jk label the incident plaquettes. Since each plaquette is shared by four edges, the factors (— l)-?i+J'2+-5'3+j4 
cancel. Thus the edge normalizations become (— l)^YAi = l/(2i + 1). Since this is positive, we can 
multiply each vertex amplitude by v'2i + 1 to take this into account. The theta networks from (36) 
contribute a sign of (—1)^^=1 , where the ik label the edges incident on the vertex. Since each edge is 
shared by two vortices, the vertex signs also cancel. The final answer is that the vortex amplitude (36) 
becomes a product of five Wigner-Racah 6j symbols multiplied by a product of six factors of the form 

We observe that in this form the tetrahedral recoupling is equivalent to the dual amplitude formula 
first proposed by Anishetty et al. [2, 3] and later used by Diakonov and Petrov [13]. The same formula 
was used in the computational work of Dass [16, 17, 20]. It should be emphasized that previous 
derivations of this formula did not make use of the spin foam formalism. As such, the identification 
of extra labels (those not coming from original plaquettes) with intertwiners was not explicit. We 
found this distinction between plaquette and intertwiner labels to be a crucial one in constructing our 
algorithm. 

A. 3. Compatible intertwiner moves. As mentioned in Definition 2.4, when a (single) spin foam cube 

move is applied, the 12 intertwiner labels of the selected cube have to be adjusted to ensure the result 
is an admissible spin foam. In general, admissible intertwiner changes depend on plaquette and edge 
labels of the cube. However, we will show next that this dependence takes on a simple form determined 
by the edge splittings. 

For this discussion, we consider the case where the edges are split according to the ladder recoupling 
of the 18j symbol (see Figure 4), as the translation invariance of this recoupling simplifies the analysis. 
From the point of view of a cube move, edge splittings come in two kinds, see Figure 6(a). In one 
case, octahedral edges corresponding to adjacent cube plaquettes remain adjacent after the intertwiner 
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edge is introduced (type A). In the other, octahedral edges corresponding to adjacent cube plaquettes 
become bridged by the new intertwiner edge (type B). The pattern of type A and type B edges for a 
cube is shown in Figure 6(b). 

The main constraint on the admissibility of intertwiner changes is parity, cf. Definition 2.2. For 
type A splittings, the two adjacent plaquottc spins change by ±i, forcing the corresponding intertwiner 
label to change by an integer. For type B splittings, the intertwiner, now bridging the two adjacent 
plaquettes, must change by a half-integer instead. A minimal set of admissible intertwiner changes now 
consists of 1, or —1, for a type A edge, and ^ or —i, for a type B edge. Note that negative and 
positive changes should be proposed with equal probability, to satisfy detailed balance. 

In summary a spin foam cube move consists of the following: 

(1) Choose a cube in the lattice. 

(2) For each plaquette in the cube, adjust the spin independently by +i or — i. 

(3) For each type A edge in the cube, adjust the spin independently by 1, or —1. 

(4) For each type B edge in the cube, adjust the spin independently by +^ or —5. 
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